%%%%%%%%%%%%%%%2CPM功率普特征估计法%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
function r= r_snr(x)

clc;
clear all;
close all;

id=sqrt(-1);
pi=3.14159;
h=1/2;                                       % 调制指数  
M=2;                                         % 进制数  
bit_number=96*2;
symbol_number=bit_number/(log2(M));          % 码元数目
sample_number=8;                             %每码元样点个数
Tb=1/32000;
fs=1/Tb*sample_number;                       % 基带取样率
x=10;
%----------------------------------------------------------------------------------------------------------------------
% data 产生0，1比特流
data=zeros(1,bit_number);
for i=1:bit_number
    temp=rand;
    if (temp<1/2)
        data(i)=0;
    elseif (temp<1) 
        data(i)=1;
    end
end;

%----------------------------------------------------------------------------------------------------------------------
%----------------------------------------------------------------------------------------------------------------------
% Ik={+1,-1,...+(M-1),-(M-1)} M=8  将比特流映射成M进制码元
Ik=zeros(1,symbol_number);
Uk=zeros(1,symbol_number);
for i=1:symbol_number
    if (data(i)==1)
        Ik(i)=-1;
    elseif (data(i)==0)
        Ik(i)=1;
    end;
end;   
Uk=(Ik+(M-1))/2;  
%----------------------------------------------------------------------------------------------------------------------


%----------------------------------------------------------------------------------------------------------------------
% CPM g(t)=1/2LT*(1-cos(2pi*t/LT)--LRC 

L=1;
g=zeros(1,sample_number+1);
for i=0:L*sample_number
    t=i/fs;
    g(i+1)=1/(2*L*Tb);      %*(1-cos(2*pi*t/(L*Tb)));
end;
% g(sample_number+1) = 0;

q=zeros(1,L*sample_number+1);
q(1) = g(1)/fs;
for i=1:L*sample_number
    q(i+1)=q(i)+g(i+1)/fs;
end;

%----------------------------------------------------------------------------------------------------------------------
% 计算载波相位  
% 模2*pi的绝对相位值（弧度）
theta=zeros(1,symbol_number*sample_number);                                           % 绝对相位
I_signal=zeros(1,symbol_number*sample_number); 
Q_signal=zeros(1,symbol_number*sample_number);

%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
sigma_a=zeros(1,symbol_number);
fi=zeros(1,symbol_number*sample_number);
for n=1:symbol_number
   if n==1
      sigma_a(n)=0;  
   else
      sigma_a(n)=sigma_a(n-1)+Ik(n-L); 
   end
      theta_n=sigma_a(n)*h*pi;
      delta=zeros(1,sample_number);
      for m=1:sample_number
          delta(m)=2*pi*h*Ik(n)*q(m);
          fi(sample_number*(n-1)+m)=delta(m)+theta_n;
      end
end

%%%%%%%%%%%%%%% 产生加性高斯白噪声AWGN %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
snr=x;                          %dB                 
sigma=sqrt(sample_number/(2*log2(M)*10.^(snr/10)));
I_noise=zeros(1,symbol_number*sample_number);
Q_noise=zeros(1,symbol_number*sample_number);
for i=1:2:symbol_number*sample_number
    [I_noise(i) I_noise(i+1) ]=gngauss(sigma);
    [Q_noise(i) Q_noise(i+1) ]=gngauss(sigma);
end;
s_noise=I_signal+I_noise+(Q_signal+Q_noise)*id;
%%%%%%%%%%%% CPM调制信号fft变换 %%%%%%%%%%%%%%%%
I_signal=cos(fi);              % 基带I路信号
Q_signal=sin(fi);              % 基带Q路信号    
s=(I_signal+I_noise)+(Q_signal+Q_noise)*id;
s0 = (I_signal)+(Q_signal)*id;
% s=I_signal+Q_signal*id;
s_fft=abs(fft(s0));       %symbol_number*sample_number点FFT
 %figure(1);
 %plot(s_fft);
 %plot((0:length(s_fft)-1)'*fs/length(s_fft),s_fft);
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
power=zeros(1,symbol_number*sample_number);
for i=1:symbol_number*sample_number
    power(i)=s_fft(i)^2/(symbol_number*sample_number); 
    plot((0:length(power)-1)'*fs/length(power),power);
    
end
%%%%%%%计算SNR%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
P_1=0;    %接收信号0-25k频带内的功率
P_2=0;    %接收信号25-50k频带内的功率
for i=1:(25/256*symbol_number*sample_number)
    P_1=P_1+power(i); 
end 
for i=(25/256*symbol_number*sample_number+1):(50/256*symbol_number*sample_number)
    P_2=P_2+power(i);  
end
%r=abs((P_1-P_2+2*7.6923)/(P_2-2*7.6923)) ; %7.6923为25-50k频带内纯信号功率的估计值
r1=(P_1)/(P_2)*(P_1)/(P_2) 
r=10*log10(r1)
    

